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Abstract 

We reconcile two scaling laws that have been proposed in the literature for the slip length 
associated with a moving contact line in diffuse interface models, by demonstrating each to apply 
in a different regime of the ratio of the microscopic interfacial width l and the macroscopic diffusive 
length Id = (Mr/) 1 / 2 , where y is the fluid viscosity and M the mobility governing intermolecular 
diffusion. For small Id/ l we find a diffuse interface regime in which the slip length scales as 
£ ~ (IdI) 1 / 2 - For larger Id/1 > 1 we find a sharp interface regime in which the slip length depends 
only on the diffusive length, £ ~ Id ~ (Mr/) 1 / 2 , and therefore only on the macroscopic variables y 
and M, independent of the microscopic interfacial width l. We also give evidence that modifying 
the microscopic interfacial terms in the model’s free energy functional appears to affect the value of 
the slip length only in the diffuse interface regime, consistent with the slip length depending only on 
macroscopic variables in the sharp interface regime. Finally, we demonstrate the dependence of the 
dynamic contact angle on the capillary number to be in excellent agreement with the theoretical 
prediction of [J], provided we allow the slip length to be rescaled by a dimensionless prefactor. This 
prefactor appears to converge to unity in the sharp interface limit, but is smaller in the diffuse 
interface limit. The excellent agreement of results obtained using three independent numerical 
methods, across several decades of the relevant dimensionless variables, demonstrates our findings 
to be free of numerical artifacts. 
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I. INTRODUCTION 


Understanding the way in which a contact line (the line along which an interface sepa¬ 
rating two immiscible fluids intersects a solid boundary wall) moves under an imposed flow 
is a problem of fundamental importance in wetting dynamics. It underpins a wide range of 
phenomena in nature, for example in the functional adaptation of many biological systems 


-4J, as well as in technological applications: from oil recovery [5|] to inkjet printing and 
microfluidics jg]. From these examples, it has become clear that experimental progress needs 
to be accompanied by the development of physical models and computational methods that 
can accurately capture wetting dynamics, even in complex geometries. 

Any purely macroscopic treatment - in which the two fluids are considered to be separated 
by a strictly sharp interface, and subject to a strictly no-slip condition on the fluid-velocity 
at the solid wall - has long been recognised to predict a non-integrable stress singularity for 
any non-zero velocity of the contact line 7), 8]. This so-called moving contact line singularity 
is clearly at odds with the commonplace experience that the contact line does, in fact, move. 

To regularise this singularity, the macroscopic picture just described must be supple¬ 
mented by additional physics that intervenes on shorter, microscopic lcngthscales. One 
approach is to relax the no-slip boundary condition on the fluid velocity, by inserting a slip 
boundary condition in a region of microscopic size £ in the vicinity of the contact line. This 


can be done by hand 9|, [10], for example by imposing a slip velocity v s — U exp(—|:c|/£), 


where x is the distance along the wall away from the contact line, or by choosing v s = £ a/r /, 
where a and ij are the wall stress and fluid viscosity respectively. In any computational 
study, slip can also arise indirectly as a numerical artifact [ill ]. 

Another approach is to keep the no-slip boundary condition and instead to recognise 
that any interface between two fluids can never be perfectly sharp, but in practice must 
have some non-zero microscopic diffuse width l. Accordingly, the use of diffusive inter¬ 
face models jl2, 13] has gained increasing popularity in recent years. A diffuse interface 
model of liquid-gas coexistence then accommodates contact line motion by a condensation- 
evaporation mechanism which, in some small region within close proximity of the contact 
line, slowly transfers matter from one side of the interface to the other 14MlT]. Similarly, 
a diffuse interface model of coexisting immiscible binary fluids accommodates contact line 
motion via a slow intermolecular diffusion of the two fluids across the interface between 
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them, again acting in a small region in the vicinity of the contact line 


Diffuse interface models 


12 


B \m, 


18|. 


| ll3] are also convenient in computational practice. An order 
parameter (“phase field”) is introduced to distinguish one fluid from another, with the length 
scale that characterises the variation of this field across the interface between the fluids pre¬ 
scribing the interfacial width, l. An advection-diffusion equation (the Cahn-Hilliard equation 
supplemented by an advective term) then determines the evolution of this order parameter, 
and hence of the interface position. This is solved in tandem with the Navier-Stokes equa¬ 
tion for the fluid velocity field. A key advantage of this approach is that all computational 
grid points are treated on an equal footing, without any need to explicitly track the position 
of the interface. In this way, even complex flow geometries can be simulated conveniently. 


It is important to emphasize, as first elucidated by jlj], that although the physical origins 
of contact line motion may differ according to the detailed microscopic physics invoked (wall 
slip, intermolecular diffusion across the interface, etc.), the macroscopic hydrodynamics far 
from the contact line nonetheless converges to a universal solution that is informed by the 
microscopies only in the sense of being rescaled by a single parameter, the “slip length” £, 
that emerges from this underlying microscopic physics. An important problem within any 
model of the microscopies is therefore to determine this emergent slip length £: not only 
because it is a key variable that determines the macroscopic wetting and fluid dynamics, 
but also because it sets a scale to which other length scales (surface heterogeneities, droplet 
size, etc.) must be compared. 

With that motivation in mind, the primary aim of this work is to determine the scaling of 
the slip length within a diffuse interface model of immiscible binary fluids. We assume the 
fluids to have matched viscosity 77 , and denote by M the mobility parameter characterising 
the rate of intermolecular diffusion in the Cahn-Hilliard equation. I 11 the existing literature, 

B least) three apparently contradictory scalings have been proposed. Several authors (e.g. 

, Bl) suggest that £ ~ Id where Id = (M 77) 1 / 2 is the characteristic lengthscale below 
which intermolecular diffusion dominates advection and above which the opposite is true. 
For convenience we call this the diffusion length in what follows, and we emphasise that 
it is determined only by the macroscopic quantities M and 77 . In contrast, [18] suggest a 
different scaling, £ ~ a /IdI, which is the geometric mean of the macroscopic diffusion length 
and the microscopic interface width. And while these two proposed scalings do not have 
any dependence on the imposed flow velocity, | 20 ] in contrast suggested that the slip length 
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does depend on the velocity as £ oc V]£ l,/2 . 

In view of these apparently contradictory predictions, an important contribution of this 
work is to provide a coherent picture that fully reconciles two of these existing predictions, by 
showing each to apply in a different regime of parameter space, and that carefully appraises 
the third. This is achieved by extensive numerical studies performed across imprecedently 
wide ranges of the two relevant control parameters. Specifically, we explore four decades 
of the dimensionless ratio Mr) /l 2 = I‘d/1 2 of the (squared) macroscopic diffusion length to 
the (squared) microscopic interface width; and three decades of the capillary number Ca 
(defined below), which adimensionalises the imposed flow velocity Vo in terms of an intrinsic 
interfacial velocity scale. Moreover we use three different, independently coded numerical 
methods, and show their results to be in excellent quantitative agreement across all decades. 

For small values of the capillary number Ca, corresponding to low imposed flow velocities, 
we demonstrate that the slip length converges to a well-defined value that is independent 
of the velocity V 0 . Any dependence on the velocity is only observed for Ca > 0.01, and we 
comment on the extent to which our results are consistent with the prediction £ oc V 0 7 of 
20] in this fast flowing regime. 

The remainder of the paper focuses on the slow flow regime, Ca < 0.01, in which the 
slip length £ is independent of the flow rate Vo- Here we show that, in the limit of large 
Id/ l, the slip length scales as £ oc Id, informed only by the macroscopic lengthscale l D - This 


agrees with the original prediction of 


13j, [ll| and corresponds to the sharp interface limit 


of the diffuse interface model. In this limit the microscopic length l essentially drops out of 
the problem, apart from providing a “singular perturbation” to regularise the contact line 
singularity. In contrast, for small Id/ l we find that the slip length scales as £ oc (Z^Z) 1 / 2 , as 


proposed by jn|. This corresponds to the diffuse interface limit in which the emergent slip 
length does depend on the underlying microscopic length l. In this way we reconcile these 
two previously apparently contradictory scalings, by showing each to apply in a different 
limiting regime of Id/ l- The crossover between the two is furthermore consistent with the 
onset of the sharp interface regime discussed by 19], 

The fact that the slip length £ is in general different from the interface width /, and indeed 
greatly exceeds it for large Id/ l, has a clear practical consequence for any simulation: the 
resolution of the full wetting dynamics appears only on a lengthscale corresponding to the 
larger of the interface width and the slip length. This, in turn, limits the region of parameter 
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space that can be considered reliable in any numerical study. 

We provide additional evidence for the distinction between the diffuse interface limit 
l D /l « 1, in which the microscopic physics informs £, and the sharp interface limit Id/1 1, 
in which it does not, by performing further simulations in which the we modify the interfacial 
(microscopic) contribution to the underlying free energy functional of the binary fluid. We 
show that the slip length depends on this modification only in the diffuse interface regime. 
In the sharp interface regime of large Id/ l, the slip length continues to depend only on 
the macroscopic dynamical variables, ij and M, free of this modification to the microscopic 
interfacial term. 

Finally, we test the validity of Cox’s formula [1] for the dependence of the dynamic contact 
angle on the Capillary number in the different slip length regimes. Our numerical results 
are in good agreement with Cox’s analytical result if we allow the slip length to be rescaled 
by a dimensionless parameter. Moreover this parameter appears, suggestively, to converge 
to unity in the sharp interface limit, but is smaller in the diffuse interface limit. 

This paper is organized as follows. In the next section we describe the models, methods 
and setups we employ to compute the slip lengths and dynamic contact angles. We present 
our results in section IIVI Finally, we summarise our key findings and discuss avenues for 
further work in section El 


II. MODEL, GEOMETRY AND BOUNDARY CONDITIONS 

In this section we specify the model that we shall use throughout the paper, starting 
with the thermodynamics in Secs. Ill Al and III Bl then the dynamical equations of motion in 
Sec. Ill Cl We then specify the flow geometry and boundary conditions in Sec. Ill Dl 


A. Thermodynamics 


We consider a binary mixture of two mutually phobic fluids, A and B , and denote the 
volume fraction of fluid A by the continuum phase held 0(r,f). The volume fraction of 
fluid B is then simply 1-0 by mass conservation and need not be considered separately. We 
consider a Landau free energy 21] 


T = 


/ dV = 

'V JV L 


r 1 m 2 

|(0 2 -l) 2 + ^f(V0) 2 


dV, 


( 1 ) 
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which allows a coexistence of two bulk phases: an A-rich phase with (f>A = 1 and a B-rich 
phase with (f>B — — 1. The bulk constant G has dimensions of energy per unit volume 
(or equivalently of force per unit area, and so modulus). If the two fluids have different 
affinities to the solid surface, a surface (wetting) contribution to the free energy can be 


added 


3,G, 


23 ] to the right hand side of Eq. (JT]). Throughout this work we assume 


neutral wetting conditions, for which no such contribution is needed. 

The chemical potential /i follows as a functional derivative of the free energy density ^ 
with respect to 0, giving 

ix = -G<1> + G(p 3 - Gl 2 V 2 0- (2) 

In equilibrium the chemical potential /x = 0. Assuming a flat interface of infinite extent in 
the y — z plane, with the interfacial normal in the x direction and the interface located at 
x = 0, we then obtain an interfacial solution of the form 

x 


= tanh 


V21 


( 3 ) 


with a homogeneous B-rich phase for x <C — l and homogeneous A-rich phase for 1 > 1. The 
interfacial constant l specifies the characteristic length scale over which 0 varies in between 
these two phases, and so corresponds to the interfacial width. The surface tension associated 
with this interface is given by 


a = 


2V2 Gl 


( 4 ) 


B. Generalisation of the Landau Theory: Introduction of a curvature term 


So far, we have specified the Landau free energy in the form most commonly used in the 
literature. It will also be instructive in what follows to consider the extent to which our 
results for the slip length do or don’t depend on the microscopic details of the model used. 
Accordingly, we now generalise the free energy slightly to give 


tf' = / vl)' h dV = 

J V JV L 


G /±2 n2 Gl 2 . \2 Gl 4 2 2 

-(0 2 - l) 2 + — (V0) 2 + a— (V 2 0) 2 


dV. 


( 5 ) 


Compared to the original free energy defined above, this has an additional interfacial curva¬ 
ture term of amplitude set by a. The mapping between this form of Landau free energy and 
the (Helfrich) continuum elastic energy can be found in e.g. 24]. The associated chemical 
potential is 


fJL — G 


+ 


l 2 V 2 0 - al 4 V 


( 6 ) 
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It is important to note that the bulk terms are unaltered, with the modification affecting 
only interfacial gradient terms involving powers of /V. A careful comparison of the original 
model, for which a = 0, with this generalised model, for which a > 0, will allow us to 
demonstrate that the slip length is independent of the microscopic details, as specified by a, 
in the sharp interface regime ljj/l 1. In contrast in the diffuse interface regime Id/ l “C 1 
we find that the slip length does depend on the microscopies, via a. 


C. Equations of Motion 


e.g. 


re dynamics of the order parameter held is specified by the Cahn-Hilliard equation [see 
2l| generalised to include an advective term: 


(dt + v.V) </> = MV 2 /n 


(7) 


Here M is the molecular mobility, which we assume constant. The fluid velocity and pressure 
fields, v(r ,t) and p(r, t), obey the continuity and Navier-Stokes equations 


d t p + v.Vp = -pV.v, 
p (d t + v.V) v = 7/V 2 v — Vp — 0V/J. 


( 8 ) 

(9) 


We denote by p and r/ the fluid density and viscosity respectively, assuming throughout that 
the two fluids are perfectly matched in both density and viscosity. In the Navier-Stokes 
equation, gradients in the chemical potential contribute an additional forcing term to the 
fluid motion, — 0V/i [25 ]. 

Diffuse interface models are widely used in the computational fluid dynamics community 
and they have been implemented using various approaches. To ensure our results are free 
from numerical artefacts, here we have used three different numerical methods: (MI) the 
lattice Boltzmann method, (Mil) spectral method, and (Mill) immersed boundary method 
to solve the equations of motion. They are detailed in Sec. IIIII below. For Mil and Mill, we 
have assumed incompressible flow and taken the inertialess limit of zero Reynolds number 
Stokes flow. This corresponds to setting to zero the terms on the left-hand-sides of Eqn. [8] 
(incompressibility) and Eq. |9] (zero inertia). The Lattice-Boltzmann method (MI) intrinsi¬ 
cally requires a small but non-zero inertia and compressibility, though our numerical results 
confirm the effects of this difference to be negligible for the problem considered here. 


7 








Vo 


FIG. 1: Schematic illustration of typical steady-states in our channel geometry. Marked are the 
key lengthscales: slip-length £, interface displacement d, channel dimensions L x ,L y , and diffuse 
interface width l. In panel a) a shear flow is imposed by moving the plates with equal and opposite 
velocities ±Vb- In panel b) a Poiseuille flow is imposed by applying a pressure drop along the length 
of the channel, sketched here in the frame in which the fluid bridge is stationary. The definitions 
for the dynamic (9^, macroscopic) and equilibrium (9 e , microscopic) contact angles are shown. 

D. Geometry, initialisation and boundary conditions 

We consider flow between infinite flat parallel plates a distance L y apart, with plate 
normals in the y direction. In the flow direction x the cell is taken to have length L x = 2 L y , 
with periodic boundary conditions on both 0 and v. All quantities are assumed invariant in 
the z direction. The phase field is initialised in an equilibrium state with a bridge of A-rich 
phase, in which 0 = +1, separated by two vertical diffuse interfaces of width E at x = L x / 4 
and 3L X /A from B-rich phases on either side, where 0 = — 1. The equilibrium contact angle 
is taken to be 9 e = 90°, corresponding to neutral wetting conditions. Throughout we define 
the position of the interface itself by the locus 0 = 0. 

The fluid is taken to be initially at rest with v = 0 everywhere. Starting from this initial 
condition, we then implement one of two common flow protocols: boundary-driven planar 
Couette flow and pressure-driven planar Poiseuille flow. In the first of these, a constant 
shear-rate 7 is applied by moving the top and bottom plates at velocities Vo = ±\^L y . This 
deforms the phase field and for small Ca = yVo/a ;$ 0.1 a steady state is reached in which 
the interface has displaced a distance d oc ±Ca at the top and bottom walls, as shown in 
Fig. HJi. The deformed bridge is then stationary, with the contact lines moving at a velocity 
=FVo relative to the top and bottom walls. 

In the Poiseuille protocol the flow is driven by an imposed pressure drop A P along the 
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length of the channel. A steady state then develops in which the bridge migrates along 
the channel at a constant speed, with the contact lines moving along each wall with a 
measured velocity Vo- In the reference frame of the contact line, therefore, both walls move 
with a velocity —Vo- By fitting the interface in the central region of the channel (between 
y = | —> |) to a circle, see Fig. [TJo, we define the dynamic contact angle 6 d as the angle 
of intersection this circle makes with the wall. The dynamic contact angle is in general 
different from the equilibrium contact angle. 

At the plates we assume boundary conditions of no-slip and no-permeation for the fluid 
velocity. For the phase held, we implement d y /j,\ y= o t L y = 0 in method Mill. In method Mil 
it is instead more convenient to set dy<j)\ y=0t L y = d%<f>\y=o,L y — dy<f>\y=o,L y = 0 (where the 
last equality need only be imposed in the case q^O). Note that while this condition auto¬ 
matically ensures zero gradient of the chemical potential, it is actually a stronger condition 
than that in demanding all the relevant odd derivatives of 0 to vanish separately. We do, 
however, find no difference in our numerical results between these. For method MI, the 
bounce-back boundary conditions assure we have no slip and no permeation of either fluid 
across the boundary. 

E. Definition of the slip length 

Our definition of the slip length uses the fact that the slip mechanism in the binary fluid 
model is intermolecular diffusion. In steady state, the Cahn-Hilliard equation of motion for 
the order parameter reads 

v.V0 = MV 2 /i. (10) 

In Fig. [21 we plot the way in which v.V0 typically varies along the fluid-fluid interface 
(taken as the locus 0 = 0, as noted above). We take the slip length £ to be given by the 
distance from the wall to the local maximum. In measuring this typical distance away from 
the wall over which v.V0 = MV 2 /j remains appreciable, £ characterises the lengthscale over 
which intermolecular diffusion is important. 

Other definitions of the slip lengths have been used in the literature to characterize the 
diffusive mechanism. For examples, [ji| defined the slip length as the point at which MV 2 // 

reaches -10% of the value at the wall after passing the first maximum, as measured along the 

□ 

interface; [19| suggested the use of the location of the stagnation point in the flow. As we 
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FIG. 2: Illustration of our definition of the slip length £, defined as the position of the first local 
maximum of v.V</> as measured along the interface (defined as the locus <fi = 0). 


will show in the results sectioiijWe recover the same scaling laws for the diffuse and sharp 


interface limits as in 


18] and [19] respectively, strongly suggesting these definitions are not 


independent, as expected if there is a genuine characteristic length scale in the problem. 


F. Parameters and dimensionless groups 

For the model and flow geometry that we have now defined, we have the following pa¬ 
rameters: the fluid density p, the fluid viscosity rj, the modulus G that sets the scale of the 
free energy density, the mobility M, the diffuse interfacial width l, the box dimensions L x 
and L y , and the applied flow velocity V 0 . (Also relevant numerically are the spatial mesh 
size and the timestep; we have ensured convergence for all the results presented.) As noted 
above, the fluid density is set to zero (inertialess flow) for methods Mil and Mill, and small 
for MI. This leaves p, G, M, /, L y , L x and Vo as the relevant parameters. We are then free to 
choose units of mass, length and time, thereby reducing by three the number of parameters 
that must be explored numerically. Accordingly, we work in what follows with the following 
dimensionless groups: y/Mrj/l, l/L y , L x lL y and the capillary number Ca = r)V 0 /a. The 
last of these characterises the speed of the imposed flow compared to the intrinsic interfacial 
velocity scale a/r] formed by comparing the surface tension to the fluid viscosity (recall that 
the surface tension a oc Gl). Among these we set the aspect ratio of the box L x = 2 L y 
throughout. The ratio l/L y is always set small, to ensure that the results are not contami¬ 
nated by finite box size effects. This leaves y/Mrj/l and Ca as the two dimensionless groups 
to be explored numerically. 
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III. NUMERICAL METHODS 


In this section we describe onr three numerical methods. For convenience we will refer to 
these throughout the manuscript as MI, Mil and Mill. We emphasize that each was coded 
independently of the others, one by each of the authors of the manuscript, thereby providing 
confidence that our results are free of algorithmic artefacts. 


A. Method I: Lattice Boltzmann Method 


We use a standard free energy lattice Boltzmann method to solve the binary fluid 
equations of motion. The basic idea behind the lattice Boltzmann algorithm is to as¬ 
sociate distribution functions, f(r, t) = {/,(r, t)} and g(r,t) = {g^r.t)}, discrete in 
time and space, to a set of velocity directions e*. Here we use a D2Q9 model, where 
ej = (0, 0), (±1, 0), (0, ±1), (±1, ±1). The physical variables are related to the distribution 
functions by 

P ^ ^ fii P^a ^ ^ fi&ion 4 * ^ ^ 9ii 4 ^ a ^ ^ 9i^ior (H) 

i i i i 

The time evolution equations for the distribution functions can be broken down into two 
steps. In the collision step, we have used a single relaxation time approach for the order 
parameter and a multiple relaxation time approach for the fluid density, 


f'(r, t) = f(r, t) - M _1 SM[f(r, t) - f eq (r,t )}, 

g'(r, t) = g(r, t) - i [g(r, t) - g e? (r, t)}. (12) 


As shown by 


23 ], the multiple relaxation time approach is significantly more reliable to 


reduce spurious velocities at the contact line. f eq and g e<? are local equilibrium distribution 
functions. The choice of the local equilibrium functions must recover the correct thermo¬ 
dynamics and hydrodynamics equations of motion in the continuum limit. The matrix M 
performs a change of basis to more physically relevant variables, including the density, mo¬ 
mentum, and viscous stress tensor. The matrix S is a diagonal matrix and contains the 
information about how fast each of these physical variables relaxes to equilibrium at every 


time step. Detailed expressions for f eg , g eq , M, and S can be found in 


23|. 
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The relaxation parameter r for the viscous stress tensor in S can be related to the 
kinematic shear viscosity v = 77 /p via 


v= (At(r — l/2))/3 . 


(13) 


Similarly, the relaxation parameter r^, for the order parameter is related to the mobility 
parameter in the Cahn-Hilliard equation through 

M = Air A - if . (14) 

T is a constant and At is the simulation time step. In practice, the relaxation times can 
therefore be tuned to match the desired continuum dynamical variables, M and 77 . 

In the propagation step, the updated distribution functions are passed on to the neigh¬ 
bouring lattice points, 


f(r + e,Af, t + At) = f'(r, t ), 

g(r + ejAf, t + At) = g'(r,t). (15) 


At the two walls we use a bounce-back rule for both the {/*(r, t)} and {^(r, f)} distribution 
functions. These ensure boundary conditions of no slip, no permeation, and no diffusion for 
the fluid material across the boundary 26]. 


B. Method II: Spectral Method 


Within this method, at each numerical timestep we separately (a) solve the hydrodynamic 
sector of the dynamics to update the fluid velocity field v at fixed phase field </>, then (b) 
update the phase field at fixed fluid velocity field. In part (a) we use a streamfunction 
formulation to ensure that the incompressibility condition is automatically satisfied. After 
eliminating the pressure by taking the curl of the generalised Stokes equation (Eqn. [9] with 
the left hand side set to zero) we solve for the streamfunction using a Fourier spectral 
method. 


For the phase field dynamics in (b) we use a third order upwind scheme [27] to update the 
convective term on the left hand side of Eqn. [3 The gradient terms on the right hand side 
are solved using a Fourier spectral method using both sine and cosine modes in the periodic x 
direction and only cosine modes in the y direction for consistency with the imposed boundary 
conditions d y (j) = d y (j) = d y (f> = 0 at the walls y = 0, L y . 
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For the time-stepping we adopt a method originally proposed by [28(, and later studied 


in depth by 29]. This splits the free energy term into an expansive part, —V 2 0, and a 
contractive part, V 2 (0 3 ), —V 4 0. These are then distributed between the n th and (n + l) th 
timesteps as 


Ln+l 


At 


= MG 


(-( 2 v>+ 2 vy)” +1 + (vy 3 - 3 vyy' 


(16) 


where At is the timestep. This method permits significantly larger timesteps than, for 
example, a Crank-Nicolson algorithm jso| in which all terms are split equally between the 
n th and (n + l) th timesteps. When present, for a > 0, the sixth order gradient term is 
treated fully implicitly. 


C. Method III: Immersed Boundary Method 


In methods I and II, the walls of the flow cell were included by means of a simulation 
box that is closed in the flow gradient direction y. In method III we instead use a biperiodic 
box that conveniently enables us to solve both the hydrodynamic sector of the dynamics 
(again in the incompressible stre am fun ction formulation at zero Reynolds number) and the 
diffusive part of the concentration dynamics using fast Fourier transforms in both spatial 
dimensions. 

To incorporate the walls of the flow cell at y — 0, L y , we then include a set of immersed 
boundary forces using smoothed (Peskin) delta functions as source terms in the Stokes 
equation along the desired location of each wall, as discussed in [3l| and [32]. The force 
required to ensure zero velocity at the location of each delta function is then evolved as 

dtF = —k (V — V 0 ), (17) 


where Vo is the prescribed wall velocity. Here k is a spring constant, and all the results 
shown below have been converged to the limit k —> oo in successive runs. 

For the phase field, we assume the boundary conditions d y /j = 0 at the plates y = 0, L y . 
This is implemented in an analogous way to the no-slip boundary condition, by means of 
adding an extra source term contribution 


d t (f) = —A"n.V/i (18) 

to the equation of motion for the concentration at the location of the Peskin delta functions, 
where n is a unit vector normal to the wall. 
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IV. RESULTS 


We now present our results. We start in Sec. IIV Al by considering the standard Landau 
free energy with a = 0, before in Sec. IIV ESI appraising the robustness of our findings to the 
inclusion of the additional interfacial gradient (curvature) term in the free energy, setting 
a > 0. In Sec IIV Cl we compare our numerical results for the dynamic contact angles to the 
analytical predictions by [l|. 


A. Standard Landau theory 


As discussed previously, in the existing literature two apparently contradictory scaling 
aws have been proposed for the contact line’s slip length £. Several authors, including 
13] and [19!'], have proposed that the slip length is proportional to the diffusive lengthscale 
Id, which describes the length below which intermolecular diffusion dominates advection 
and above which the opposite is true, giving £ oc In — (Mr/) 1 / 2 . In contrast, the lattice 
Boltzmann simulations and scaling arguments of [l8| suggest that the slip length depends 
not only on this macroscopic diffusive lengthscale Id, but also on the microscopic interfacial 
width /, with £ oc (Z^/) 1 / 2 . 

To resolve this apparent discrepancy we have performed extensive numerical simulations 
of the moving contact line problem across four decades of the relevant dimensionless control 
parameter Mr//l 2 = l 2 D /l 2 - Moreover, to ensure that our results are free of algorithmic 
artefacts we have used (across the entire range) the three different, independently designed 
and coded numerical techniques just described. We focus the discussion in this section on 
the case of planar Couette flow, returning to consider planar Poiseuille flow later in the 
manuscript. We also focus on the slow flow limit Ca —> 0, returning later on to consider the 
effects of finite Ca. 

The results are shown in Fig. [3] As can be seen, the three numerical techniques give 
virtually indistinguishable results across all four decades of Mr]/1 2 . Over these four decades 
we can distinguish two distinct regimes: a (i) diffuse interface regime when Mr]/l 2 <C 1, in 
which £ oc (/jV) 1 / 2 , an d (h) a sharp interface regime when Mr]/l 2 1, in which £ oc Id- In 
between these two distinct regimes is a broad crossover window that itself spans around the 
decade either side of Mr]/l 2 ss 0.2. In this way, importantly, our results encompass and unify 
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FIG. 3: Data for the slip length £ from all three numerical methods for planar Couette flow clearly 
demonstrating both the diffuse- and sharp-interface limits for small and large Mrj/l 2 respectively. 
All the results here are in the slow flow regime, Ca < 0.01, where the slip length is independent of 
the shear velocity. 


both of the two previously apparently contradictory scalings put forward in the literature, by 
showing each to apply in a different regime of the relevant dimensionless control parameter 
Mr]/l 2 . 

In addition to the two basic lengthscales Id and l present in the model equations, out 
of which the slip lengthscale £ emerges in the manner just described, there are two other 
lengthscales present in our simulations, set by the system size L y and the discretisation scale 
Ax. We have ensured that the results presented in Fig. [3] are independent of possible finite 
size effects due to L y and discretisation errors due to Ax. The former is a particular hazard 
in the limit of large Mr]/l 2 , where the slip length £ becomes large, and the latter in the limit 
of small Mr]/l 2 , where the slip length becomes small. 

Having shown our numerical results to capture both the sharp interface limit of £ oc Id 
for large Mr]/l 2 , and the diffusive interface limit of £ oc (IdI) 1 ^ 2 for small Mr]/l 2 , we are now 
in a position to reprise carefully the scaling arguments put forward in the earlier literature 
for each of these two scaling forms separately, and to discuss the validity of the assumptions 
made in arriving at these functional forms. 

We will follow a similar line of arguments as in jl9|, where the Stokes and steady-state 
Cahn-Hilliard equations are integrated in the x-direction along the wall. For concreteness, 
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let us focus on the contact line at the top wall, i.e. at y — L y . The same scaling analysis 
applies for the bottom wall (y — 0). For the ^-component of the Stokes equation, this gives 


d 2 v x d 2 v x 


dp 

dx 


dp 

dx 


dx = 0. 


(19) 


dy 2 dx 2 

The second term can be integrated to r\dv x jdx^^ = 0. The contribution of the third term 
in the integral is also zero, because the pressure attains the same, constant value on either 
side of the interface, far from the contact line. The fourth (final) term scales as /i max , which 
represents the magnitude of the chemical potential for a given flow setup. 

The first term in Eq. [T9] needs to be analysed more carefully. In the sharp interface limit, 
the term d 2 v x /dy 2 scales as Vo/£ 2 , and this term is significant across a length scale given by 
the slip length £. Thus, we have 


/^max ' 


' Mmax? 

vVo/S- 


( 20 ) 


In the diffuse interface limit, the fluid velocity still varies across a length scale £ in the 
y-direction, and as such, the term d 2 v x /dy 2 still scales as Vo/£ 2 . However, to capture the 
key physics, the spatial window of integration must be broadened to the interface width l, 
since here £ is less than the interfacial width l : recall the bottom leftmost data points in 
Fig. [3] Taking the suitable of window of integration into account, Eq. [T9] leads to 


h £2 l ~ pmajx.i 
/Xmax ~ hW/£ 2 ‘ 

Let us now carry out a similar analysis for the Cahn-Hilliard equation, 


dx Vy dy 


M 


9V 


M- 


dx = 0. 


( 21 ) 


( 22 ) 


dx 2 dy 2 

The first term can be integrated to give Vo d'|^° 00 = 2Vo- The second term is zero simply 
because v y = 0 at y = 0. The third term integrates to —Mdp/dx 1^, which is again zero 
because dp/dx = 0 far from the interface. 

The result for the fourth term depends on the slip length regimes. d 2 p/dy 2 scales as 
P max/£ 2 - However, the window of integration is different in the sharp (£) and diffuse (/) 
interface limits. For the sharp interface limit, the Cahn-Hilliard equation leads to the scaling 


Vo ~ Tf/i max /£. 


(23) 
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On the other hand, for the diffuse interface limit, we obtain 


fd) ~ M/w^ 2 - (24) 

We are now in the position to derive the scaling law predictions for the slip length. 
Combining Eqs. (liZUji and (T2WI) gives 

f ~ (Mr]) 1 / 2 = l D (25) 

in the sharp interface regime. Similarly, combining Eqs. (}2Tj) and (I24p for the diffuse interface 
regime results in 

f ~ (Mrf) 1 /^ = (l D l) 1/2 . (26) 

As seen in Fig. El in between the sharp and diffuse interface regimes is a smooth crossover 
where the behaviour varies smoothly from one scaling law to another. We define the crossover 
point by fitting the power laws in the sharp and diffuse interface regimes separately, and 
finding where these two fits intersect each other. This occurs at I n/I = 0.38, in broad 
agreement with the sharp interface criterion l^/l > 0.25 proposed by 19] 

We note finally that the values of the slip length presented in Fig. 0 are independent 
of the wall velocity, and correspondingly of the Capillary number Ca = 77 V 0 /CT, provided 
Ca < 0.01. (Data not shown.) We shall demonstrate this independence explicitly below for 
the case of planar Poiseuillc flow (see Fig. [5]). 

B. Influence of a curvature term 

In this section we test the extent to which our results for the slip length do or do not 
depend on the microscopic details of the diffuse interface model. Intuitively we might expect 
the slip length to be independent of microscopies in the sharp interface regime, but dependent 
on microscopies in the diffuse interface regime. Our numerical results below provide evidence 
in favour of this intuition, to which we shall however also add a note of caution at the end 
of this section. 

There are many possible ways to modify the basic Landau theory for which we presented 
results in the previous subsection. We focus here on the simple but non-trivial extension 
made by introducing the additional interfacial curvature term with a strength set by a in 
Eq. (J5]). Fig. E] shows our numerical results for a = 0.1 and 1.0. Shown in the same plot 
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Mr/fl 2 


FIG. 4: Numerical data demonstrating how slip lengths in the diffuse-interface limit are affected by 
the inclusion of a curvature free energy term. The results of the standard Landau theory (a = 0) 
are shown for comparison. Inset: percentage change relative to the case without curvature for 
a = 0.1 (red squares), a = 1 (green circles). These data are obtained using the spectral method 
(Mil) in planar Couette flow. 

for comparison are our original results, already presented in the previous subsection, for the 
original theory with a = 0. As can be seen, the introduction of the curvature term affects 
the slip length quite strongly in the diffuse interface regime. This is to be expected: in 
this regime the physics of the problem is determined not only by macroscopic quantities, 
but by the microscopic gradient contributions to the free energy. In contrast, as the control 
parameter Mr)/l 2 increases into the sharp interface regime the dependence of the slip length 
on this microscopic parameter a dramatically decreases and appears to become negligible 
in the sharp interface limit Mr)/l 2 —>■ oo. 

Although these numerical results strongly suggest that the slip length is independent of 
the microscopic (gradient) contributions to the free energy in the sharp interface regime, 
we attach the following note of caution. In the case of equilibrium phase coexistence in 
the absence of flow, it is possible to show by integrating across the interface between the 
phases that macroscopic bulk quantities, such as the value of the chemical potential at phase 
coexistence, are independent of the structure of that interface, as determined by the spatial 
gradient terms in the free energy. However the same is not true in general for the case 
of phase coexistence under the non-equilibrium conditions of an applied flow because the 
equations of motion are usually non-integrable in that case: they cannot be integrated across 
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the interface to give a result that depends only on the bulk quantities on either side of the 
interface, but the result instead retains a dependence on the structure of the interface itself. 
This issue has been discussed in particular detail in the case of stress selection in shear 


banded flows 


33 


34]. Therefore, while our numerical results do strongly suggest that the 


slip length is independent of the gradient terms in the sharp interface limit, this question 
merits some further attention in future studies. 


C. Relation to Cox’s result 

The results presented so far have demonstrated that, in the sharp interface limit, the 
emergent slip length is determined by the macroscopic dynamical variables 77 and M. While 
this is clearly an important finding in terms of our physical understanding of the moving 
contact line problem, in numerical practice this sharp interface limit Mr]/l 2 —> 00 can be 
difficult to attain. Indeed the time taken to attain a steady state after the inception of flow 
increases dramatically with increasing Mr]/l 2 . Lengthy simulations were required to obtain 
the rightmost data points in Fig. [3l with the extreme cases requiring run times of the order 
of weeks with a parallelised code using 4 cores of Intel(R) Xeon(R) 2.40Ghz. 

With this in mind, we now turn to address two questions of practical numerical impor¬ 
tance: (i) whether simulations carried out in the diffuse interface regime still reproduce the 
expected macroscopic dynamics far from the contact line, and (ii) whether our definition of 
the slip length is still meaningful outside the sharp interface regime. We shall address these 
questions in the context of the Poiseuille flow protocol sketched in Fig. [1 }d, in particular 
asking whether the dependence of the observed dynamic contact angle on flow velocity is 
the same in the diffuse interface regime as in the sharp interface regime. 

We start by checking that the slip length obtained for the Poiseuille flow protocol is 
equivalent to that for planar Couette flow. Fig. Or, shows that the slip length for the 
Poiseuille case converges to a constant value as Ca —>■ 0, as for the planar Couette case. We 
have checked in this regime that the measured values of the slip length for the Couette and 
Poiseuille flows are within 0.5% of each other over the whole range of Mr]/l 2 . See the inset 
of Fig. [5}i, where the slip lengths for the Poiseuille flow are compared to the best fit power 
laws extracted for planar Couette flow in both the diffuse and sharp interface regimes. For 
larger Ca > 0.01 the slip length starts to depend on the fluid velocity, as in planar Couette 
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FIG. 5: Numerical data for the Poiseuille flow protocol, for Mg/l 2 = 0.01 (black circles 0)> 0-1 (red 
squares), 1 (green diamonds), 4 (blue triangles), a) Demonstration that the slip length converges 
to a constant value as Ca —>• 0. The inset shows that the measured value of using the 

Poiseuille flow protocol is in good agreement with the best fit power laws obtained for shear flow 
in Fig [3l b) Symbols mark the dynamic contact angle 6 d as measured from simulation. The solid 
lines show the analytical predictions of [1] where 6 = (3£(j a /L y . The inset shows the best fit values 
of /3 which appear to converge /3 — > 1 in the sharp interface limit, Mg/l 2 — >• oo. The simulation 
data are obtained using the spectral method (Mil). 


flow. 

We now turn to the dynamic contact angle, as defined in Fig. [TJo. It is known that this 
angle increases with increasing flow speed: for increasing pressure drop A P in Fig. [T}d the 
interface is deformed more significantly due to the larger flow-speed mid channel compared 
to the stationary wall, and correspondingly the dynamic contact angle deviates more from 
the equilibrium value. Associated with this is a larger gradient in the chemical potential, 
which indeed allows the contact line to move at a higher velocity. 

An analytical prediction for this increase of the dynamic contact angle with flow speed, 
as characterised by the capillary number Ca , was given by {l|: 

g(9 d )=g(9 e ) + Ca\n(5- 1 ). (27) 

Here 9 d is the dynamic contact angle and 9 e is the contact angle at equilibrium (here 90°). 
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For fluids of matched viscosity the function g is defined as 


9(0) = 


7T0(7T — 0) + (2,77 0 — 7T 2 ) sill 0 COS 0 — 7T sill 2 0 
2n 2 sin 0 


d(j). 


(28) 


The parameter 5 in Eqn. [23 is the ratio of the microscopic slip length to some characteristic 
macroscopic lengthscale of the system. We choose that to be the channel width L y and set 
5 = /3^(j a /L y where 0 is a constant fitting parameter. 

To test Cox’s prediction, we plot our simulation results for the dependence 9d on the 
capillary number Ca in Fig. [5b. We do so for four different values of Mr)/l 2 , which we 
ensured span both the diffuse and sharp interface limits, as well as the crossover regime 
between them. Pleasingly we find excellent agreement with Cox’s prediction across the the 
full range of Ca, for all values of Mr)/l 2 . Furthermore, as shown in the inset of Fig. [5b, the 
fitting parameter 0 appears suggestively to converge to unity in the sharp interface limit, 
Mr)/l 2 —>■ oo. 

These results suggest that simulations in the diffuse interface limit (and similarly the 
crossover regime) still provide a reliable representation of the macroscopic dynamics. How¬ 
ever, the effective slip length must be suitably corrected by an 0(1) prefactor 0. Increasingly 
minor corrections are needed as one approaches the sharp interface limit, and we expect Cox’s 
result would be directly obtained for Mr)/l 2 —>■ oo, where 0 —> 1. 


V. CONCLUSIONS 

In this paper, we have performed extensive numerical simulations to reconcile previously 
apparently contradictory scaling laws for the slip length associated with a moving contact line 
in diffuse interface models. We did so by exploring four decades of the relevant dimensionless 
control parameter Mg/l 2 = l 2 D /l 2 , and three decades of the capillary number Ca. Moreover 
we used three independent numerical methods, demonstrating their predictions to be in 
excellent quantitative agreement across all decades. 

In the limit of slow flows corresponding to typical capillary numbers Ca < 0.01, we 
find that the slip length converges to a well-defined value that is independent of the flow 
speed. In this slow flow regime we explored four decades in Mg/l 2 = l 2 D /l 2 . Doing so 
enabled us to distinguish two distinct regimes for the scaling of the slip length, and a 
broad crossover window in between. For Mr)/l 2 C 1 we found a diffuse interface regime 
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in which the slip length scales as £ ~ (IdI ) 1//2 , consistent with the previous result of 18], 
and is furthermore sensitive to form of the microscopic interfacial terms in the free energy 
functional. In contrast for Mr]/l 2 > 1 we recover the sharp interface limit previously 


discussed by [13] and [191 ]. Here the slip length is proportional to the diffusion length scale, 
£ oc l D ~ (M??) 1 / 2 , depending only on the macroscopic dynamical variables of viscosity r] and 
molecular mobility M. Our numerical results further suggest the slip length to be insensitive 
to changes in the microscopic interfacial free energy terms in this sharp interface regime. 

For Ca > 0.01, the effective slip len gth decreases with increasing capillary number, in 


qualitative agreement with the work of [20]. However, we are not able to reproduce their 
prediction £ oc V 0 7 of in this fast flowing regime. More broadly, so far we are imable to 
find any reliable scaling law that is valid across a broad range of Capillary number, and for 
both the sharp and diffuse interface regimes. 

Our numerical results also allow us to appraise suitable sets of simulation parameters. By 
comparing our simulation results for the dynamic contact angle to the analytical prediction 
by 1|, we find excellent agreement if we allow the slip length to be rescaled by a dimensionless 
prefactor (3 = 0(1), which pleasingly appears to converge to unity in the sharp interface 
limit. In numerical practice, however, the sharp interface limit can be very expensive to 
realise, partly because the time taken to attain a steady state increases dramatically with 
increasing Mr]/l 2 , and partly because the slip length is 0(10) larger than the interface 
width, with demanding implications on the overall size of the simulation box. In the diffuse 
interface limit, in contrast, the timescale and lengthscale requirements are less stringent, and 
simulations accordingly much less expensive to perform. Our fits to the Cox formula suggest 
that a reliable representation of the macroscopic physics is nonetheless attainable even in 
the diffuse regime, provided the lengths are suitable rescaled by this corrective prefactor 
(3 = Oil). 

There are several avenues for future work. Firstly, it would be interesting to investigate 
how the two scaling regimes for contact line slip in diffuse interface models affect more 


complex flow problems, such as the macroscopic dynamics of a liquid droplet 
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criterion for the liquid entrainment (wetting failure) [38|, [39], and whether similar renormal¬ 


37 | or the 


ization of the slip length is adequate to account for their differences, as for the Poiseuillc 
flow considered here. In these cases, we need to implement non-neutral wetting boundary 
conditions. We expect the prefactors of the scaling laws to change with contact angles, but 


22 














not the exponents themselves. 

Secondly, given the crossover we have reported in this study, it would be insightful to 
revisit the liquid-gas system and understand the sharp interface limit of this model. In 
the binary model the slip mechanism is diffusion, and as such, the slip length depends on 
the mobility parameter, in addition to the fluid viscosity and the interface width. In the 
liquid gas system, the mechanism is evaporation-condensation instead, and we expect the 
slip length to depend on the density ratio between the liquid and gas phases. 

Thirdly, we note that we have taken the zero temperature limit in this work and have not 
accounted for thermal fluctuations. It will be interesting to study if such fluctuations simply 
broaden the effective interface width, or they give rise to new phenomena. For example, 


enhanced spreading due to thermal fluctuations has been reported by 


40] and 41]. 


Finally, experimental systems consisting of colloid-polymer mixtures are known to phase 
separate into colloid-rich and colloid-poor domains, with an interface width that is of order 


1 micron 




much larger than the typical value for molecular fluids. Such a system 


can potentially be exploited to realise the different contact line slip regimes discussed here 
experimentally. 
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